i'JAv. 

MOirr:s.vt.. 












I'l* 



f I 

«> 





NAVAL POSTGRADUATE SCHOOL 



Monterey, California 



APPROXIMATE IM ERVAE ESTIMATION 
METHODS EOR HIE 

REEIABII ITYOE SYSTEMS I SIXG COMPONENT 

DA I'A WITH exponp:ntial 

AND WEIBl'LL DISTRIBUTIONS 








by 



LEE, I lyeon-Soo 



September, 1989 



Thesis Advisor 



W.M. WOODS 



Approved for public release; distribution is unlimited. 



nclassified 

:uriiy classification of this pacie 



RKPORI DOCUMLMAIION RAGi: 



1 Report Security Classification L.'nclassified 


lb Restrictive Markings 


Security Classification Authority 


3 Distribution Availability of Report 

Approved for public release; distribution is unlimited. 


3 Declassification Do\\ngradinp Schedule 


Performing Organization Report Numberfs) 


5 Monitoring Organization Report Number(s) 


1 .Name of Performing Organization 
.'aval Postgraduate School 


ob Office Symbol 
(if applicable) 55 


7a Name of Monitoring Organization 
Naval Postgraduate School 


: Address (clrv, siaie, and ZIP code) 

lonierey. CA 93943-5000 


7b Address (cir\\ siaie, and ZIP code) 

Monterey. CA 93943-5000 


1 Name of Funding Sponsoring Organization 


8b Office Symbol 
( if applicable) 


9 Procurement Instrument Identification .Number 


: Address (clry\ siaie. and ZIP code) 


10 Source of Funding Numbers 


Program Element No Project No Task No W ork Unit Accession .No 



; Title (include securir.-ciassijicaiion) APPROXIMATE INTERVAL ESTIMATION METHODS rOR THE RELIABILITY 
SYSTEMS USING COMPONENT DAI A WITH EXPONENTIAL AND WEIBUl L DISTRIBUTIONS 



l Personal Author(s^ EEE. Hvcon-Soo 



la Type of Report 


13b Time Covered 


14 Date of Report (year, month, day ) 




15 Page Count 


laster's 1 hesis 


From To 


September, 19S9 




64 



) Suppleinentary Notation I'he view s expressed in this thesis are those of the author and do not reflect the official policy or po- 
tion of the Department of Defense or the U.S. Govermnent. 



' Cosati Codes 


18 Subject Terms [continue on lexerse if necessar\ and Identify by block number) 


leld 


Group 


Subgroup 


Reliability, Jacknife, Confidence Limit, ExponentiakWeibuU.Nomal, NewTon - Raphson 



) Abstract \conilnue on r€\ersc ij necessary and ideniify by block ?rumber) 

Two approximate parametric inter\'al estimation methods for system reliability using component test data are developed 
id evaluated. One method can be applied to any coherent system with components which have exponential failure limes 
ith possibly dilTerent failure rates and different mission operating times. 1 his method estimates the ratios of component 
lilure rates which are then used to de\ clop the approximate lower confidence limit. These ratio estimates are developed with 
id without jacknife methods and the two results are compared. This procedure is veiy accurate and simple to compute, 
‘quiring the use of standard chi-square tables. 1 his ratio method is subsequently extended to coherent systems with com- 
onents whose failure times have a Weibu// distribution. A nearly exact parametric lower confidence limit for P{X > x) is 
sveloped and evaluated where a is given and A’ has a normal distribution with unknown mean and variance. Ehis procedure 
also simple to evaluate and requires the use of Student l tables. 



) Distribution .A\ailability of .Abstract 
3 unclassified unlimiied □ same as report 
Name of Responsible Individual 

XM WOODS 

D FORM 14"3,sa MAR 



□ DTIC users 


21 Abstract Security Classification 
L’nclassified 




22b Telephone ( include Area code) 

(408)454-1893 


22c Office Symbol 

55Wo 


83 APR edition may be used until exhausted 
All other editions are obsolete 


security classification of this page 






L'nclassified 



1 



1245264 



Approved for public release; distribution is unlimited. 



Appro.ximate Interval Estimation .Methods for the Reliability 
of Systems Using Component Data with Exponential and 
WeibuII Distributions 

by 

LEE, 1 Iveon-Soo 

Captain, Republic Of Korea Army 
B.S., Korea .Military Academy, 1982 

Submitted in partial fulfillment of the 
requirements for the degree of 

MASfER OF SCIENCE IN OPEICATIONS RESEARCH 

from the 

NAVAL POSIGICADUATE SCHOOL 
September. 1989 



ABSTRACT 



Two approximate parametric interval estimation methods for system reliability using 
component test data are developed and evaluated. One method can be applied to any 
coherent system with components which have exponential failure times with possibly 
different failure rates and different mission operating times. This method estimates the 
ratios of component failure rates which are then used to develop the approximate lower 
confidence limit. These ratio estimates are developed with and without jacknife methods 
and the two results are compared. This procedure is ver\’ accurate and simple to com- 
pute, requiring the use of standard chi-square tables. This ratio method is subsequently 
extended to coherent systems with components whose failure times have a U’cibull dis- 
tribution. A nearly exact parametric lower confidence limit for P{X > ,v) is developed and 
evaluated where ,v is given and A' has a normal distribution with unknown mean and 
variance. This procedure is also simple to evaluate and requires the use of Student t ta- 
bles. 
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I. BACKGROUND 



Parametric confidence interval procedures for the reliability of mechanical systems 
are much less developed than procedures for electrical systems. This is due to the more 
complicated failure distributions used to model mechanical hardware. The failure rate 
of a series system of independent components, with exponentially distributed failure 
times and equal mission times, is the sum of the failure component failure rates. This 
property has permitted development of numerous methods for system reliability of series 
and other coherent systems using component failure data. The IVeibull and extreme 
value distributions have been used in life testing methods for both electrical and me- 
chanical devices. Several methods have been used for obtaining point estimates of pa- 
rameters for these distributions and for the reliability function itself, Harter and Moore 
[Ref 1: pp. 889-901], Mann [Ref 2; pp. 231-256]. 

The derivation of simple confidence limits for the reliability function for the extreme 
value distribution with parameters c and 6 has posed problems, because methods based 
on a pivotal quantity such as [z — z)jd are inadequate. Lawless [Ref 3; pp. 355-364], 
Johns and Lieberman [Ref 4: pp. 135-175] and Thoman. Bain and Antle [Ref 5: pp. 
363-372] have developed nearly exact procedures for confidence limits for the reliability 
function of the U'eibull and extreme value distributions. Schneider and Weissfeld [Ref 
6: pp. 179-186] have developed interval estimation methods for percentiles of the U'eibull 
and extreme value distribution based on censored data. Although extensive methods 
have been developed for interval estimates of the reliability of a single component with 
the U'eibull or extreme value failure distribution, very few parametric interval methods 
have been developed for system reliability using component test data with U'eibull failure 
distributions. 

Two approximate parametric interval estimation methods for reliability of coherent 
systems using component test data are developed and evaluated in this thesis. Evalu- 
ations of these procedures are performed using computer simulation for series systems 
only. One of these methods is developed for the reliability of a series system whose 
components have exponential failure distributions and difl'erent mission times. This 
procedure was found to be quite accurate and can be applied to coherent systems in 
general. This first procedure is then extended to the case where components of the sys- 



tern have Weibidl failure distributions. The method used is an extension of a non- 
parametric method developed by Myhre, J., Rosenfeld, A. and Saunders, S. [Ref 7; 
pp.21?-227j. 

The normal distribution is used extensively in some mechanical rcliabilty models. 
Maximum likelihood and minimum variance unbiased estimators for P{X>x^) when 
both and a- are unknown were developed over thirty years ago, Lieberman and 
Resnikoff [Ref 8: pp. 457-516], Folks and others [Ref 9: pp. 43-50]. and Barton [Ref 
10: pp. 227-229]. Exact interval estimation procedures for P{X>x^) were developed by 
Owen and Hua [Ref 11: pp.285-311] using the non-central t distribution. Letting X de- 
note strength and Y denote stress, mechanical reliability is sometimes modeled as the 
value for P{X > )'). Approximate interval estimation procedures for this probability 
when A' and F are assumed to be normal have been developed by Church and Harris 
[Ref 12: pp. 49-54] when the mean and variance of }' are known. Downton [Ref 13: pp. 
551-558] modifies their procedure slightly to get more accurate bounds and suggests an 
approximate procedure when the means and variances of both X and }' are unknown. 
Lower confidence intervals for P{X> Fj obtained under the assumption of normality for 
X and }' can lead to serious error when either X or }' or both are truncated well into the 
tails. Consequently, P(X > .\y) may be a more reasonable model of mechanical reliability 
where .v„ is chosen conservatively. A very accurate approximate lower confidence limit 
procedure for P(A'>jrJ is developed and evaluated in this thesis. It can be computed 
easily. 
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II. INTERVAL ESTIMATION PROCEDURE - EXPONENTIAL CASE 



A system of independent components is coherent if an increase in reliability of any 
one of its components does not cause a degredation in system reliability. Suppose a 
coherent system has k components. We assume that the failure distribution of compo- 
nent ; is exponential with failure rate . Then system reliability can be written as a 
function of / = 1, 2, • • • . A; i.e., 



^ 5(0 ~ ^- 2 ’ ' " ^ ^-k' hi ^ h<) 

where /, = /,(/) is the operating time for component /, / = 1. 2. • • • . . 



Let /.„ = max{/i. /j, • • • , /,*} and r, = / J/.„ / = 1, 2. • • • . A. Then one can write 

^s(0 = sUmi hi'" 1 *'ki hi • • * < h<) 

If the were known and were an upper 10()(1 — a.y/o confidence limit for , the 
corresponding lower confidence limit for R^{t) would be 

r A 

h- • • • i>'k- h-'" 1 he) 

Specifically, if the system is a series system of independent components, so that 



then. 




( 2 . 1 ) 



If n, items of component / are tested until failure, T denotes the total test time accumu- 

A A 

lated by all n, items and n = V/;,, then is y\„ . See Bain and Engelhardt (Ref 

1=1 1=1 

14 ). 

An upper confidence limit for is 
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( 2 . 2 ) 



2 

2n 

^■/n,C'(x) ~ j. 




i=] 



where yX„ is the 100(1 — a)th percentile of the yl distribution. Corresponding equations 
for truncated testing are similar. 

If the r, are unknown, the following methods estimate the values for r, from the data. 
One method uses the likelihood ratio estimate for r,. The second method uses a jacknifed 
version of the first method. The two resulting confidence limits Rs,u^) Rs,lm . 
and without jacknifing r, respectively, are compared for relative accuracy. Quenouille 
[Ref 15; pp. 353-360] first reported a method for estimating ratios that reduced the bias 
without increasing the variance. Miller [Ref 16: pp. 1-15] gives an excellent review of the 
jacknife method which includes a discussion on the application of jacknifing to estimat- 
ing ratios. 

A. LOWER CONFIDENCE LIMIT WITHOUT JACKNIFING 

In this case, the maximum likelihood estimate of the ratio /• = is 



where a, = hJT, and = max (/i,. •••,/*) . The resulting approximate upper confi- 
dence limit for is 



A 



k 

where n = Tn, , and 
/=! 




X 



2 



2n 



k 







= 1 . 2 , 



,k 



(2.4) 



(2.5) 



and T, denotes the failure time of the ^th test for component /. 
The resulting appro.ximate confidence bound Rs,u,) is given by 
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( 2 . 6 ) 



R 



S.L{^) ~ 




B. LOWER CONFIDENCE LIMIT Rs,, WITH .lACKNIFING 

The definitions for , r, , n, , and 7', in Section A are also used in this section. Let 

A 

denote the estimate for /I, by removing T,j from the data ; i.e.. 



and 




J = 1,2. 



(2.7) 



n„ J I, -. • • • , I J 

Vr 

/ /m/ 
t=\ 



(2.S) 



Then the jacknifed ratio estimate r' is given by 



A* » A 

>'i = « >'i - 




where «’ = min(A7,. /2,. • • • , and 



A _ '‘■iM 



j = 1,2, 



Now define 



■ C'(a) 



by 



2 

^-y,2n 

''■L\^)= ~k 




/=! 



(2.9) 



( 2 . 10 ) 



( 2 . 11 ) 



The corresponding confidence bound Rs,u-.) given by 
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( 2 . 12 ) 



R 







When the //,, i = 1,2, dilTer considerably, this jacknife estimation proce- 

dure can be unbalanced. That is, the number of data points used to compute the jacknife 
estimate will differ from one component to another. It was decided to use this jacknifing 
procedure and determine the effect of differing sample sizes by examining the results of 
the simulations. 

In equation (2.9), the jacknife estimate is constructed by using only the first n' ob- 
servations from each component to obtain where n' = min{«,. Wj, • • • , «*}. Of 
course it is rather arbitrary to take the first n' observations from , since any n' of the 
collection of n, values could be used. 



C. SIMULATIONS AND RESULTS FOR EXPONENTIAL CASE 
1. Simulation 

a. Simulation language and package 

The programming language used to simulate this problem was VS 
FORI R.'\N on an IBM 3033. Also LEXPN in LLIUANII was used to generate observed 
exponential random variates. SHSORT was used to perform the sort routines. 

b. Cases and Input parameters 

The six input parameters below determine the conditions for exercising the 
simulation runs. System reliability is determined by five parameters. Values of these 
parameters are gi^ en in the tables that show the results. The test plan simulated was to 
test n, items until all fail. 

number of component types, k : 5 and 15 
system reliability. : 0.9 and 0.975 

significance level, a : 0.05 and 0.2 

component time, /, : small to large (see tables) 

sample size, ii, : small to large (see tables) 

k 

Reliability of a series system is expressed as Rs = exp( — jV . We chose arbitrarily 
to determine the failure rates. A, , from this equation by assuming all to be equal, 
consequently. 
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- In Rs 
A/, 



c. Replications 

The procedure was replicated 1000 times for each case to get 1000 values 
o\' Rs,nx) snd 1000 values of Rs,n^) ■ ^ ^ order each set of the 1000 values of Rs.u^) ^nd 
^ 5 ,Lw iri ascending order. Then the two approximate confidence bounds. Rsxwimm--^)) 

A 

(1000(1 ->)) ' ihe 1000(1 — a)th order values of these two sets of data. Finally, 

for each of the two ordered sets of data, we find the order indices y, and for which 

A 

f^s.Lwu\) ^s.Lwu 2 ) closest to Rs for their respective sets. Then y, / 1000 and 

J 2 1 1000 are called the two corresponding simulated true confidence levels. 

2. Results 

Tables 1 through 3 show the results of the 3 cases simulated. The results indi- 
cate that this interval estimation method using estimates of failure rate ratios will yield 
quite accurate lower confidence limits for system reliability when components have un- 
known constant failure rates. The jacknife method also yields very accurate confidence 
limits which are slightly conservative. That is the 100(1 — a) percentile points of Rs,lm < 
given in the tables, are slightly less than the true value of Rs- Consquently, 
— 1 — a . Alternatively, one can say Rs.lm ^ conservative 100(1 — a) 

percent lower confidence limit procedure for Rs . 



7 



Table 1. RELIABILITY OF A SERIES SYSTEM WITH SMALL NUMBER 



( LESS THAN 10 ) OF SAMPLE SIZES - EXPONENTIAL CASE 



Number of 
Compo- 
nents 


Reliability 


Cf. 


Lower Confidence 
Limit 


True Confidence 
Limit 


of System 


W O 

Jacknifing 


WITH 

Jacknifing 


WO 

Jacknifing 


WITH 

Jacknifing 




.90 


2 


0.9017 


0.8904 


0.7720 


0.9100 


5 


.05 


0.8985 


0.8889 


0.9620 


0.9780 


.975 


• ^ 


0.9750 


0.9720 


0.8020 


0.9160 




.05 


0.9746 


0.9729 


0.9600 


0.9800 




.90 


. X. 


0.9019 


0.8901 


0.7600 


0.9500 


1 s 


.05 


0.9012 


0.8933 


0.9300 


0.9800 


1 0 


.975 




0.9745 


0.9724 


0.8300 


0.9600 




.05 


0.9746 


0.9725 


0.9600 


0.9800 



Sample sizes for 5 components are 9. 7. 10, 8. 6 

Sample sizes for 15 components are 6. 7. 5. 6. 7. 8. 9. 5. 8. 7. 9. 10. 7. 9. 6 
Component times i, are 2. 3. 7. 8. 10, 3. 7. 10, 1, 7, 8, 3, 10, 1, 8 



Table 2. RELIABILITY OF A SERIES SYSTEM WITH MEDIUM NUMBER 
( LESS THAN 30 ) OF SAMPLE SIZES - EXPONENTIAL CASE 



Number of 
Compo- 
nents 


Reliability 


a 


Lower Confidence 
Limit 


True Confidence 
Limit 


of System 


W O 

Jacknifing 


wnii 

Jacknifing 


W O 

Jacknifing 


WITH 

Jacknifing 




.90 


.2 


0.8999 


0.8978 


0.8020 


0.8660 


. 


.05 


0.8990 


0.8983 


0.9600 


0.9680 




.975 


2 


0.9750 


0.9745 


0.8000 


0.8460 




.05 


0.9748 


0.9745 


0.9620 


0.9740 




.90 


2 


0.8992 


0.8960 


0.8300 


0.9100 


15 


.05 


0.8980 


0.8981 


0.96O0 


0.9700 


.975 


• ^ 


0.9751 


0.9738 


0.7900 


0.9000 




.05 


0.9750 


0.9741 


0.95O0 


0.9800 



Sample sizes for 5 components are 30. 21, 10, 15. 26 

Sample sizes for 15 components are 3, 7, 10, 15, 20, 15, 7, 5, 30. 20, 10, 7, 13, 21, 30 
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Table 3. RELIABILITY OF A SERIES SYSTEM WITH LARGE NUMBER 



( LESS THAN 100 ) OF SAMPLE SIZES - EXPONENTIAL CASE 



Number of 
Compo- 
nents 


Reliabilitv 


a 


Lower Confidence 
Limit 


True Confidence 
Linut 


of System 


W 0 

Jacknifing 


WITH 

Jacknifing 


WO 

Jacknifing 


WITH 

Jacknifing 




.90 


2 


0.8999 


0.8991 


0.8080 


0.8200 


5 


.05 


0.8990 


0.8992 


0.9680 


0.9620 


.975 


.2 


0.9750 


0.9750 


0.8540 


0.8060 




.05 


0.9751 


0.9750 


0.9380 


0.9520 




.90 


2 


0.9000 


0.8988 


0.8010 


0.8 7.^0 


15 


.05 


0.9000 


0.8996 


0.9500 


0.9590 


.975 


1 


0.9 7.^0 


0.9747 


0.8460 


0.8730 




.05 


0.9750 


0.9750 


0.9520 


0.9470 



Sample sizes for 5 components are 30. 63. 75. 98. 26 

Sample sizes for 15 components are 15. 40. 35. 17, 26. 67. 50. 65, 80. 32, 95, 100, 15, 
45, 30 
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III. INTERVAL ESTIMATION PROCEDURE • WEIBULL CASE 



A. LOWER CONFIDENCE LIMIT 

Consider a scries system with k components. Let the time to failure A' of compo- 
nent / have a IVeibull distribution with density 

;:(/,) = exp{ r/'} , /,>0 . (3.1) 

Then 

;?,(/,) = exp( , /,>0 (3.2) 



and 



Rsii) = expj - ^ ;.f' /f I = expj - ;.;Y r, /f'j , / > 0 . (3.3) 

where /,' = = max/.,' and r, = If tlie /?, are known, A'f- has constant 

i 

failure rate /f- and the procedures in Chapter II can be used to obtain Rsxw and 
with T,, replaced by T^- in equation (2.5). It is well known that )'=lnA', has 
an ex i rente value distribution with CDF 

Fy{v) = I - exp| - e I 

where f=ln(l//.,) and <), = 1//?,. 

Engelhardt and Bain [Ref 17: p. 323] have developed the following simple unbiased 
estimators for c, and U, . using ordered values >7 = In 



s n 




where 5 = [0.84/?,] = largest integer < 0.84 n, and A'^,, is the yth order statistics from the 
sample of size //, of A', . Also 
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^i = ln( — )=>v + yt5 
y. 

where >’ = 0.5772 and P, = ^>y«, . Let 




(3.5) 



(3.6) 



and 



r.. = V-' 

hj 






y=i,2. 



(3.7) 



We approximate the distribution of T/j by the exponential distribution with failure rate 
/f' = and proceed as in Chapter II. Define 



= 



Ti 



where T, = i = 1. 2. •••, k. Let = max and 



. - I 



Then an approximate upper confidence limit for is given by 






. 2 

2 /? 



w,L’(a) 



k 

(=1 



(3.S) 



(3.9) 



(3.10) 



and the corresponding approximate lower confidence limit Rs,lm given by 






exp 



J 



— / 



m.Ui'x) 




(3.11) 



where n = Y^n, . 

i- 1 

This procedure is labeled the Formula procedure in the tables that follow in this 
section. Its distinguishing feature is the equation for /?, given by equation (3.6). An al- 
ternative procedure, labeled the Newton - Raphson procedure, estimates /?, using the 



II 



maximum likelihood procedure which is obtained using a Newton - Raphson approxi- 

A 

niation method. Equations for /?, under the Newton - Raphson procedure are provided 
in Appendix A. 



B. SIMULATIONS AND RESULTS 



1. Simulation 

a. Simulation language and package 

The programming language used to simulate this problem was VS 
FORTR.AN on an IBM 3033. Also LEXPN in LLRANIl was used to generate observed 
exponential random variates. SHSORT was used to perform the sort routines. 

b. Cases and Input parameters 

The six input parameters below determine the conditions for exercising the 
simulation runs. System reliability is determined by five parameters. Values of these 
parameters are given in the tables that show the results. The test plan simulated was to 
test n, items until all fail. 



number of component types, k : 5 and 15 
system reliability, : 0.9 and 0.975 

signilicance level, a : 0.05 and 0.2 

component time. /, : small to large (see tables) 

sample size, n, : small to large (see tables) 



k 



Reliability of a series system is expressed as = exp{ — V } . Failure rates /, can 

1=1 

be determined from that equation by assuming all /,/, to be equal, consequently, 






— In R^ 

~T7T 



c. Replications 

The procedure was replicated 1000 times for each case to get 1000 values 
of and 1000 values of Rs,lm • We order each set of the 1000 values of Rs,lm 

Rs,uy) in ascending order. Then the two approximate confidence bounds, /?j,a,)(iooo(i-.» ^nd 
( 1000 ( 1 -,)) , for ore the 1000(1 — (z)th order values of these two sets of data. Finally, 
for each of the two ordered sets of data, we find the order indices y, and for which 
^ 5 .z.(,)(;, ) ond Rs,LM(j 2 ) ore the closest to Rs for their respective sets. Then y, / 1000 and 
J 2 1 10*X) are called the two corresponding true confidence levels. 
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2. Results 

Tables 4 and 5 display the results of the simulations and determine the accuracy of 
ill equation (3.1 1) as a lower confidence limit procedure for system reliability, 
Rs. for parameter values /„ /!„ r, and Rs given in the tables. The terms in the tables have 
the same meaning as the corresponding terms in Tables 1 through 3 which were e.x- 
plained in Section 2.C.2 . The procedure would be exact for the Formula method if the 
values in the Formula column equal the corresponding numbers in the same row in the 
Reliability of System column. In Table 4 for example, the 80th percentile point of .9205 
for the Newton - Raphson procedure is more accurate than the Formula procedure 
which has an 80th percentile point of .9324. The last column shows, however, that what 
we have called an 80% low^er confidence limit procedure is in fact closer to a 61% pro- 
cedure. The accuracy improves if the sample size is increased from 15 to 30 as indicated 
in Table 5. The accuracy improves even more if the number of components in the system 
increases from 5 to 15 as indicated in Tables 4 and 5. 

The results in Tables 4 and 5 show that the R^^^ method given by equation 
(3.11) is too optimistic - especially for small sample sizes and for systems with a small 
number of components. Modifications to this procedure are needed. These tables also 
indicate that the Newton - Raphson method is superior to the Formula method. 
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Table 4. ACCURACY OF R,, AS A LOWER CONFIDENCE LIMIT FOR Rs 
WITH SAMPLE SIZES OF 15 - WEIBUf L CASE 



Number of 
Compo- 
nents 


Reliability 
of System 


a 


Lower Confidence 
Limit 


True Confidence 
Limit 


Formula 


Newton 

Raphson 


Formula 


Newton 

Raphson 


5 


.90 


.2 


0.9324 


0.9205 


0.4820 


0.6060 


.05 


0.9512 


0.9406 


0.5380 


0.6660 


.975 


.2 


0.9849 


0.9824 


0.5040 


0.6540 


.05 


0.9905 


0.9891 


0.5600 


0.6660 


15 


.90 


.2 


0.9193 


0.9030 


0.5440 


0.7620 


.05 


0.9416 


0.9288 


0.6060 


0.7760 


.975 


2 


0.9802 


0.9747 


0.6160 


0.8080 


.05 


0.9865 


0.9828 


0.6580 


0.8240 



Initial beta for 5 components are 2.2. 2.4. 2.6. 2.8. 2.5 

Initial beta for 15 components are 2.2. 2.4, 2.6. 2.8. 2.5. 2.3, 2.7. 2.5, 2.2. 2.6. 2.9, 2.8, 
2.4. 2.6. 2.1 

Component times t, are 2, 3, 7, 8, 10. 3, 7, 10, 1, 7, 8, 3, 10, 1,8 



Table 5. ACCURACY OF R^, AS A LOWER CONFIDENCE LIMIT FOR Rs 
WITH SAMPLE SIZES OF 30 - WEIBULL CASE 



Number of 
Compo- 
nents 


Reliability 


a 


Lower Confidence 
Limit 


True Confidence 
Limit 


of System 


Formula 


Newton 

Raphson 


F ormula 


Newton 

Raphson 




.90 




0.9249 


0.9125 


0.4100 


0.6320 




.05 


0.9424 


0.9321 


0.5100 


0.6940 




.975 


1 


0.9837 


0.9798 


0.4380 


0.6160 




.05 


0.9879 


0.9845 


0.4600 


0.6760 




.90 


2 


0.9207 


0.9036 


0.5560 


0.7200 


15 


.05 


0.9333 


0.9198 


0.4380 


0.7520 


.975 


2 


0.9808 


0.9758 


0.5120 


0.7520 




.05 


0.9850 


0.9810 


0.5200 


0.7920 



Initial beta for 5 components are 2.2. 2.4, 2.6. 2.8. 2.5 

Initial beta for 15 components are 2.2. 2.4, 2.6, 2.8, 2.5, 2.3, 2.7, 2.5, 2.2, 2.6, 2.9, 2.8, 
2.4. 2.6. 2.1 
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IV. INTERVAL ESTIMATION PROCEDURE - NORMAL CASE 



A. BACKGROUND 

In recent years it has become popular to model mechanical reliability as P{X > h) where 
X denotes strength and Y denotes stress. This formulation of reliability is an average of 
P(.V>y| L = j') , since 



P{X> = Ey P{X > y) = 



P{X>y^fy{y^dy . 



An alternative model is the worst case approach. In the worst case model, one se- 
lects a worst case value ofy, say . and then designs the strength. A' ,of the component 
so that P(A’> yo) = Rq where Rq is a reliability requirement. If A' has a normal distrib- 
ution then this requirement imposes constraints on the mean and variance of A’ . This 
is usually done in a manner to comply with standard "safety factor" procedures. 

The average model, P{X > I"), uses two random variables and is subject to more 
random error than the worst case model. The accuracy of the expression P{X > I”) for 
values of this expression close to one is highly suspect due to deviations in the tail 
probabilities of both A' and T from those assumed in the model. It is common to assume 
that both .V and Y have normal distributions. Truncated normal distributions would be 
more appropriate for many types of mechanical equipment. Harris and Soms [Ref. IS; 
pp. 650-663] discuss implications of this problem. Very significant errors in point and 
interval estimation for reliability are readily demonstrated when A' is truncated normal 
but assumed to be normal in the more simple model which specifies that R = P{X > }). 
Table 10 shows this effect when X is truncated above at /r -f 1.645cr , where Z, is the 100 
(1 — a)th percentile point of the standard normal distribution. Church and Harris [Ref. 
12; pp. 49-54] and Downton [Ref. 13: pp. 551-55S] have developed approximate confi- 
dence intervals for P{X > }") when A' is normal with unknown mean and variance and }' 
has the standard normal distribution. 
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B. EXACT AND APPROXIMATE INTERVAL ESTIMATES 



Minimum variance unbiased estimators (MVL’E) for R = P(A'>j') are well known 
when X is normal with unknown mean and known variance and also when the variance 
is unknown. In the former case, Lieberman and Resnikofl'[Ref 8] developed the result 



A / A' - i - 

R = 0 






/ /? - 1 



\' n 



which is MVUE for R where <I> is the standard normal cumulative distribution function. 
When the variance is unknown, several versions of an integral expression for P{X>y) 
have been developed by Lieberman and Resnikoir[Ref 8] , Basu [Ref 19] and Folks and 
others [Ref. 9). 

Exact lower confidence interval estimates for P{X>y) when A’ is normal with un- 
known mean and variance involve the non-central t - distribution. Owen and Hua [Ref. 
11] developed tables of the lower 90® o and 95% confidence limit values for P(X>y) 
based on the non-central i - distribution. These values are tabled for values of k in the 
range -3.0 (.2) 6.0 and sample sizes n = 2 (1) 18, 21 (3) 30. 40 (20) 100 , where 

k = (.V — y)ls and x and s are the sample mean and standard deviation. Their tables are 
essentially exact. .An approximation to their exact tabled values is given by R'l where 

Rl = , 



= k- 



1T + 



k^ 









a , ^—1 



(4.1) 



k = (.V —y)ls and is the 100(1 — a)th percentile of the i distribution with n — 1 de- 
grees of freedom. Equation (4.1) was developed in this thesis. It is an extensive ad hoc 
modification of an equation developed by Church and Harris [Ref 12]. Tables 6 and 7 
display values of Ri = <!>(>') , y , y’, R[ . and R^ — R^^ Cor k = 1. 2. 3. 4. sample sizes n 
= 10. IS. 30 and confidence levels 90% and 95% . R^ and y denote the "exact" lower 
confidence limits and corresponding 0‘‘(RJ values from Owen and Hua [Ref. 1 1]. Both 
■■/' and R2 are given by equation (4.1). The accuracy of the approximate confidence in- 
terval is quite good relative to the values for R^ given by Owen and Hua [Ref. 1 1]. 
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Table 6. APPROXIMATE ( Rl ) AND EXACT ( R, ) 90% CONFIDENCE 
LIMITS FOR P{X> 



n 




k 


RiMy) 


y 


y' 


7^L:d>(y-) 


Rl - R, 






1 


0.68 156 


0.47194 


0.45454 


0.67528 


-0.00628 


10 


1.3830 




0.89130 


1.23397 


1.20199 


0.88517 


-0.00613 




0.97453 


1.95262 


1.88991 


0.97025 


-0.00428 






4 


0.99602 


2.65302 


2.54950 


0.99437 


-0.00165 






1 


0.73037 


0.61385 


0.61133 


0.72950 


-0.00087 


18 


1.3334 


O 


0.92569 


1.44512 


1.44038 


0.92488 


-0.00081 


3 


0.98755 


2.24314 


2.23150 


0.98684 


-0.00071 






4 


0.99877 


3.02679 


3.00614 


0.99859 


-0.00018 






1 


0.75937 


0.70424 


0.70508 


0.75960 


0.00023 


30 


1.3114 


2 


0.94256 


1.57740 


1.57852 


0.94248 


-0.00008 


“S 


0.99233 


2.42404 


2.42459 


0.99205 


-0.00028 






4 


0.99945 


3.26207 


3.25926 


0.99940 


-0.00005 



Table 7. APPROXIMATE ( Rl ) AND EXACT ( R^ ) 95% CONFIDENCE 
LIMITS FORP(A'>}’o) 



n 


^.05.r-l 


k 




V 

/ 


/ 


/?[:(!)( -/) 


R'i-Rl 






1 


0.63052 


0.33311 


0.27702 


0.60912 


-0.02140 


10 


1.8331 


2 


0.85187 


1.04477 


0.94228 


0.82692 


-0.02405 




0.95565 


1.70308 


1.52864 


0.93655 


-0.01910 






4 


0.99031 


2.33813 


2.07743 


0.98075 


-0.00956 






1 


0.69504 


0.51007 


0.49292 


0.68896 


-0.00608 


18 


1.7396 


2 


O.9034S 


1.30221 


1.26991 


0.89777 


-O.00571 


3 


0.97996 


2.05344 


1.99739 


0.97674 


-0.00322 






4 


0.99735 


2.78717 


2.70338 


0.99638 


-0.00097 






1 


0.73361 


0.62369 


0.61789 


0.73167 


-0.00194 


30 


1.6991 


2 


0.92855 


1.46579 


1.45391 


0.92677 


-0.00178 


“S 


0.98855 


2.27522 


2.25448 


0.98758 


-0.00097 






4 


0.99894 


3.07140 


3.04028 


0.99873 


-0.0002 1 
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Tables S and 9 display the results of computer simulations with 1,000 replications 
to check the accuracy of the R[ method for 80% and 90% lower confidence liinits for 
P{X>y) s R for y = 3. with various values of o and n determined so that R equals 
the values shown. The procedure would be exact at the 80% level if the values in the 
column labelled a = .2 equal the corresponding values of R in the same row. The "true" 
confidence level corresponds to the index i{R) of 1,000 ordered values of R^,^ for which 
^ exaniple, in the seventh row of Table 8, R = .950 , a, = 20 . X = 10, 
^L(. 20 ),*oo ^^as .9575, and P-hao). 9 oo .95.?3. Also Rnjo)-^ .950 and Rn,]t))^f.t^ .950. 

Tables 8 and 9 indicate that the R[ procedure given by equation (4.1) is quite accurate 
at the 80° 0 level of confidence for the cases simulated. 
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Table 8. ACCURACY ANALYSIS OF Rf PROCEDURE FOR 80% AND 90% 
LOWER CONFIDENCE LIMITS FOR P(A>3) WHEN X IS 
NORMALLY DISTRIBUTED 



R 




n 


1 

o 

o 


True Confidence Level 




a = .2 


a = .1 


a = .2 


a = .1 






10 


0.9624 


0.9504 


0.7310 


0.8970 




0.5 


25 


0.9546 


0.9542 


0.7560 


0.8820 






75 


0.9517 


0.9501 


0.7800 


0.8980 






10 


0.9606 


0.9531 


0.7470 


0.8910 


.950 


1.0 


25 


0.9501 


0.9552 


0.7970 


0.8770 






75 


0.9509 


0.9528 


0.7900 


0.8770 






10 


0.9575 


0.9533 


0.7580 


0.8860 




20.0 


25 


0.955S 


0.9551 


0.7410 


0.8730 






75 


0.9511 


0.9511 


0.7830 


0.8820 






10 


0.9924 


0.9913 


1 .0000 


0.8840 




0.5 


25 


0.9906 


0.9S9S 


0.7800 


0.9070 






75 


0.9903 


0.9907 


0.9800 


0.S7S0 






10 


0.9921 


0.9S99 


0.9990 


0.9000 


.990 


1.0 


25 


0.9913 


0.9904 


0.7570 


0.8910 






75 


0.9904 


0.9S98 


0.7800 


0.9060 






10 


0.9923 


0.9922 


0.9880 


0.8750 




20.0 


25 


0.9909 


0.9908 


0.7760 


0.8 800 






75 


0.9902 


0.9903 


0.7940 


0.8920 






10 


0.9966 


0.9948 


1 .0000 


0.903O 




0.5 


25 


0.9955 


0.9949 


0.7790 


0.9020 






75 


0.9956 


0.9949 


0.7540 


0.9060 






10 


0.9966 


0.9952 


1.0000 


0.8950 


.995 


1.0 


25 


0.9960 


0.9955 


0.8630 


0.8880 






75 


0.9949 


0.9947 


0.8040 


0.9150 






10 


0.9970 


0.9951 


1.0000 


1.0000 




20.0 


25 


0.9960 


0.9957 


0.7720 


0.8700 






75 


0.9951 


0.9948 


0.7940 


0.9140 
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Table 9. ACCURACY ANALYSIS OF Rf PROCEDURE FOR 80% AND 90% 
LOWER CONFIDENCE LIMITS FOR P(Y>30) WHEN X IS 
NORMALLY DISTRIBUTED 



R 




n 


^L.iooon->) 


True Confidence Level 




a = .2 


0! = .1 


a = .2 


a = .1 






10 


0.9621 


0.9466 


0.7510 


0.9000 




0.5 


25 


0.9530 


0.9519 


0.7680 


0.8890 






75 


0.9507 


0.9491 


0.7880 


0.9070 






10 


0.9605 


0.9528 


0.7480 


0.8910 


.950 


1.0 


25 


0.9497 


0.9549 


0.8030 


0.8810 






75 


0.9509 


0.9527 


0.7910 


0.8780 






10 


0.9575 


0.9533 


0.7580 


0.8860 




20.0 


25 


0.9558 


0.9551 


0.7410 


0.8730 






75 


0.9511 


0.951 1 


0.7830 


0.8820 






10 


0.9923 


0.991 1 


1 .0000 


0.8860 




0.5 


25 


0.9900 


0.9890 


0.7990 


0.9140 






75 


0.9901 


0.9904 


0.7980 


0.8860 






10 


0.992 1 


0.9899 


0.9990 


0.9310 


.990 


I.O 


25 


0.9911 


0.9902 


0.7620 


0.8950 






75 


0.9904 


0.9S9S 


O.78O0 


0.9070 






10 


0.9923 


0.9922 


0.9880 


0.8750 




20.0 


25 


0.9909 


0.9908 


0.7760 


0.8800 






75 


0.9902 


0.9903 


0.7940 


0.8840 






lo 


0.9965 


0.9946 


1 .0000 


0.904C» 




0.5 


25 


0.9952 


0.9946 


0.7900 


0.9140 






75 


0.9955 


0.9948 


0.7620 


0.9100 






10 


0.9966 


0.9951 


1 .0000 


0.9280 


.995 


I.O 


25 


0.9959 


0.9953 


0.7610 


0.8920 






75 


0.9949 


0.9947 


0.8040 


0.9140 






10 


0.9970 


0.9951 


1.0000 


1.0000 




20.0 


25 


0.9960 


0.9957 


0.7600 


0.8700 






75 


0.9951 


0.9948 


0.7940 


0.9340 
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Table 10 displays the inaccuracies of the R[ lower confidence limit procedure for 
P(X > 3) using equation (4.1), which assumes A' is X{n. o-), when in fact A’ has the dis- 
tribution of a normal random variable with mean a and variance a- that has been trun- 
cated at fx + 1.645(7. Note that n -f 1.645(7 > > 3, because P{X> 3) = .95, .99 and .995. 

Examination of Table 10 reveals gross inaccuracies; consequently, even when the 
distribution of A' is truncated far above the value y, the exact lower confidence limit for 
/’(A'> r) can be greatly in error when computed under the assumption that A' is normal. 
This problem will be compounded when one is computing “exact" confidence intervals 
for T(A’> }") assuming both X and Y are normal when in fact one or both may be irun- 
catcd normal. This suggests that modeling mechanical reliability as P{X > }') may be 
more risky than more standard engineering approachs for modeling mechanical reliabil- 
ity which include the notion of safety factors. 
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Table 10. ACCURACY ANALYSIS OF Rf PROCEDURE FOR 80% AND 90% 
LOWER CONFIDENCE LIMITS FOR P(A > 3) WHEN X IS TRUN- 
CATED NORMAL 



Reli- 
abil- 
ity of 
Sys- 
tem 


<7, 


n 


^L,1000(l-«) 


True Confidence Level 


a — .2 


a = .1 


a = .2 


a = .1 






10 


0.9806 


0.0566 


1.0000 


1.0000 




0.5 


25 


0 . S 941 


0.0125 


1.0000 


1.0000 






75 


0.9383 


0.0871 


1.0000 


1 .0000 






10 


0.8904 


0.0419 


1.0000 


1.0000 


.950 


1.0 


25 


0.9555 


0.0813 


1 .0000 


1 . 0()00 






75 


0.9431 


0.0782 


1 .0000 


1 .0000 






10 


0.8503 


0.0105 


1 .0000 


l.OOOO 




20.0 


25 


0.9274 


0.0485 


1 .0000 


1 . 0 ( too 






75 


0.9496 


0.1222 


1 .0000 


1 .0000 






10 


0.9734 


0.6035 


1 .0000 


1.0000 




0.5 


25 


0.9924 


0.6346 


1 .0000 


1 .0000 






75 


0.9780 


0.6552 


1 .0000 


1 .0000 






10 


0.9976 


0.9172 


1 .0000 


1 .0000 


.990 


1.0 


25 


0.9884 


0.7505 


1.0000 


1 .0000 






75 


0 . 995 1 


0 . 680 O 


1 .0000 


1 .0000 






10 


0.9593 


0.5116 


1 .0000 


1 .0000 




20.0 


25 


0.9989 


0.6821 


1 .OOOO 


1 .0000 






75 


0.9883 


0.7514 


1 .0000 


1 .0000 






10 


0.9813 


0.7328 


1 .0000 


1 . 0 () 0 ( ) 




0.5 


25 


0.9719 


0.9096 


1 .0000 


I.OOOO 






75 


0.9933 


0.9130 


1.0000 


1 .0000 






10 


0.9709 


0.9015 


1.0000 


1 .0000 


.995 


1.0 


25 


0.9991 


0.8996 


1.0000 


1.0000 






75 


0.9949 


0.9267 


1.0000 


1 .0000 






10 


0.9723 


0.6432 


1 .0000 


1 .0000 




20.0 


25 


0.9923 


0.7539 


1.0000 


1.0000 






75 


0.9971 


0.9489 


1 . 00()0 


1.0000 



V. CONCLUSIONS AND RECOMMENDATIONS 



The lower interval estimation procedure for reliability of coherent systems which 
was developed in Chapter 11 appears to be accurate, easy to use, and applicable to co- 
herent systems. Although this procedure assumes that failure times of all components 
of the system are independent and have exponential probability distributions, it can be 
easily extended to systems with component failure distributions that can be transformed 
into exponential failure distributions ; e.g., WeibuU distribution with known shape pa- 
rameter. This procedure has potential for being combined with a similar procedure for 
systems with cyclical components. The combined procedure would provide for use of 
binomial component test data and exponential component test data to compute lower 
confidence liiriits on the reliability of coherent systems with both cyclic and continuous 
components. Such an extension could use a failure rate ratio estimation procedure sim- 
ilar to that de\ eloped in this thesis. Such a method should be explored. 

The interval estimation method for the reliability of a system with components that 
have WeibuU failure times is not sunwiently accurate to be applied to systems that have 
10 or fewer components each with ten or fewer tests. This procedure needs further study 
and refinement. 

The approximate lower confidence limit for component reliability P{X>y) when A’ 
is normally distributed with unknown mean and variance is very accurate. It does not 
require an extensive set of tables such as those developed by Owen and llua[Ref 11], 
but only requires the use of the standard t tables. 
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APPENDIX A. MLE OF WEIBULL PARAMETERS BY THE NEWTON - 

RAPHSON METHOD 

Let X ~ /?). Then the likelihood function for the first r ordered observations 

from a random sample of size n is given by 



• • • • Xr.n) = 



n: 






( 



i=\ 



Cl-fAfe,)]' 



= 4 Tri( ^ - [Z( ^ - 'X ^ )' 

(=1 •- i=l 



Setting the partial derivatives of this likelihood with respect to a and P equal to zero 
gives the MLE's a and P as solutions to the equations 




n 




and 



A 



a = 




/=i 



n 



1 



/? 



where n = r. It can be shown that these equations have unique solutions which are the 
maximum likelihood estimates. The NEWTON - IL.'XPHSON method for solving an 

A A 

equation gip) = 0 is to determine successive approximations p,. where 
/?„, = p, — giPPjg'iP). Therefore, the estimates of a and P can be solved by letting 
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oooooononnnonnoonoooooooononoooonoooooooooooooon 



APPENDIX B. FORTRAN CODE FOR INTERVAL ESTIMATION 
PROCEDURE - EXPONENTIAL CASE 

PROGRAM EXPONE 



V:*Vc*Vryc’Vf'5VyfVf:VVfVf*Vf*VfVoV^V‘}VVf*VfVf‘5VVrVfVf*Vf*Vc'*VfVfVf**Vf'>V’5V'}V;V'>V*3V‘5V'5VVc‘5V')V'>V‘5V'>V'>V'»V*’5VVc'5VVfVcVr'>V*3V‘5V 



* 



Vr 



* THIS IS A PROGRAM TO COMPUTE THE TRUE CONFIDENCE LIMIT OF * 

* A SERIES SYSTEM WITH / WITHOUT JACKNIFE METHOD AND * 

* COMPARE THE DIFFERENCE OF THOSE RESULTS. * 

Vf * 



* BELOWS ARE GIVEN OR ASSUMED. 



* N ; NUMBER OF COMPONENTS TYPE * 

* RSYS ; RELIABILITY OF A SERIES SYSTEM * 

* ALPHA ; SIGNIFICANCE LEVEL * 

* TIME ; TEST TIME * 

* SAMPLE ; SAMPLE SIZE FOR EACH TYPE OF COMPONENT * 

* *jV 

Vr Vf 

* THESE ARE VARIABLES USED. * 

* * 



* 



* 




B I GLAM ; 
CASE ; 
CHISQR ; 
CLRl ; 
CLR2 
CTIME ; 
DF 

DIFF ; 
LAMBDA ; 
LAMHAT ; 
LAMHM ; 
LAMHST ; 
LAMMAX ; 
LOMIT ; 
KEYl ; 
KEY2 ; 



R1 

R2 ; 

RATIO 1 ; 
RATI02 ; 
ROM IT ; 
ROMSUM ; 
RSTAR ; 
RVALl ; 
RVAL2 ; 
T ; 

ZALPHA ; 



TEMPORARY ARRAY FOR EXPONENTIAL RANDOM VARIATE * 
LARGEST VALUE OF LAMBDA 

COUNTER * 

CHI-SQUARE VALUE FOR GIVEN ALPHA AND DF 

TRUE CONFIDENCE LIMIT FROM NON- JACKNIFING * 

TRUE CONFIDENCE LIMIT FROM JACKNIFING * 

TEST TIME OF INDIVIDUAL COMPONENT 

DEGREE OF FREEDOM FOR CHI-SQUARE. * 

ABSOLUTE VALUE OF DIFFERENCE BETWEEN R1 & R2 * 

FAILURE RATE OF EACH COMPONENT TYPE 

LAMBDA HAT FOR JACKNIFING 

FINAL LAMBDA HAT FROM NON -JACKNIFING * 

FINAL LAMBDA HAT FROM JACKNIFING * 

LARGEST VALUE OF LAMHAT * 

LAMBDA WITH 1 COMPONENT OMITTED 
TEMPORARY ARRAY * 

TEMPORARY ARRAY * 

INITIALLY COMPUTED RELIABILITY FOR EACH 
COMPONENT TYPE * 

COMPUTED RELIABILITY BY NON-JACKNIFING 
COMPUTED RELIABILITY BY JACKNIFING * 

RATIO OF LAMBDA FOR NON-JACKNIFING * 

RATIO OF LAMBDA FOR JACKNIFING 

RATIO WITH 1 COMPONENT OMITTED * 

USED FOR JACKNIFE, SUM OF R WITH OMIT 1 
FINALLY COMPUTED R VALUE ( R STAR ) * 

R( 500 * (1-ALPHA) ) FOR NON-JACKNIFING * 

R( 500 * (1-ALPHA) ) FOR JACKNIFING * 

TOTAL TEST TIME OF EACH COMPONENT TYPE 
RIGHT PERCENTILE POINT( NORMAL DISTRIBUTION) * 
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VfVr*Vc**Vf***Vf*Vrx*VrVrVf ************■>>* Vf'5V*'5Vy(rVf’>V’5V’>V'5V'jV**Vc'5V*Vc’^VVc'jVVc'jV/V*Vc'>VVcVc’VcVcVfVc’3VVc’ 



PARAMETER (NN = 500) 

REAL TIME(15), ALPHA(2), ZALPHA(2), RSYS(2), LAMBDA(15), P 
REAL LAMMAX, LAMHAT(15), B( 100) , RSTAR(15), LAMHM 

REAL RATI0K15), RATI02(15), BIGLAM, DIFF(15) 

REAL L0MIT(15,100), CTIME( 15 , 100) , R0MIT( 15 , 100) , T(15) 

REAL ROMSUM , SUMl, SUM2 , LAMHST, Rl(NN), R2(NN), KEyi(NN) 
REAL RVALK15), RVAL2(15), CLR1(15), CLR2(15), KEY2(NN) 

REAL DIFFRK15), DIFFR2(15), R(15) 

INTEGER SAMPLE(15), N(2), DF, CASE 

DATA N / 5 , 0 / 

DATA RSYS / .9, .975 / 

DATA TIME / 2,3,7,8,10,3,7,10,1,7,8,3,10,1,8 / 

DATA SAMPLE / 30,63,75,98,26,15,7,5,30,20,10,7,13,21,30 / 

DATA ALPHA / .2, .05 / 

DATA ISEED / 1736 / 



/-• SUBROUTINE ZTABLE COMPUTES THE RIGHT PERCENT POINT ZALPHA */ 
/■'■ FROM RIGHT CUMULATIVE PROBABILITY ALPHA */ 

DO 10 I = 1, 2 

CALL ZTABLE (ALPHA( I), ZALPHA(I)) 

10 CONTINUE 



CASE = 0 

II IS INDEX FOR N */ 

DO 150 II = 1, 2 

COMPUTE THE DEGREE OF FREEDOM FOR CHI-SQUARE */ 



DF = 0. 

DO 20 I = 1, N(II) 

DF = DF + SAMPLE(I) 
20 CONTINUE 

DF = 2 * DF 



/* FINDING COMPONENT TYPE THAT HAS THE MINIMUM NUMBER */ 

/'■• OF SAMPLE SIZE */ 

NSTAR =999 

DO 25 I = 1, N(II) 

IF (SAMPLE(I) .LE. NSTAR) THEN 
NSTAR = SAMPLE(I) 

ENDIF 

25 CONTINUE 
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/* JJ IS INDEX FOR RSYS 
DO 140 JJ = 1, 2 



*/ 



/* COMPUTE LAMBDA FROM THE GIVEN EQUATION AND */ 

/* FIND THE BIGGEST LAMBDA */ 

BIGLAM = 0 

DO 30 K = 1, N(II) 

LAMBDA(K) = ( -AL0G(RSYS( JJ) ) / N(II)) / TIME(K) 

IF( LAMBDA(K) . GE. BIGLAM) BIGLAM = LAMBDA(K) 

R(K) = EXP( - LAMBDA(K) * TIME(K) ) 

30 CONTINUE 



/* COMPUTE CHI-SQUARE VALUE FOR GIVEN ALPHA AND DF */ 



/* LL IS INDEX FOR ALPHA */ 

DO 130 LL = 1, 2 

IF ( DF .EQ. 1 ) THEN 
P = ALPHA(LL) / 2 
CALL ZTABLE(P, ZALPHA(LL)) 

CHISQR = ZALPHA(LL) ** 2 

ELSE IF ( DF .EQ. 2 ) THEN 

CHISQR = -2 * ALOG(ALPHA(LL)) 

ELSE IF ( DF .GE. 3 ) THEN 

CHISQR = DF * (1 - 2./(9 * DF) + 

ZALPHA(LL) * SQRT( 2./(9 * DF))) ** 3 

END IF 



CASE = CASE + 1 
DIFF(CASE) = 0 

DO 120 L = 1, NN 

LAMMAX = -99. 

DO 50 I = 1, N(II) 



/* GENERATE EXPONENTIAL RANDOM */ 

/* NUMBERS WITH MU = 1 */ 

CALL LEXPNdSEED, B, SAMPLE(I), 1, 0) 

T(I) = 0 

DO 40 J = 1, SAMPLE(I) 

/* CONVERT TO EXPONENTIAL RANDOM */ 
/* NUMBERS WITH MU = LAMBDA AND */ 
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40 

50 

60 

70 

80 

90 



/* ADD THOSE FOR EACH COMPONENT */ 
/* TYPE */ 

B(J) = B(J) / LAMBDA(I) 

T(I) = T(I) + B(J) 

CTIME(I,J) = B(J) 

CONTINUE 

LAMHAT(I) = SAMPLE(I) / T(I) 



/* FINDING MAXIMUM LAMBDA HAT AND ITS */ 
/* INDEX */ 

IF ( LAMHAT(I) . GE. LAMMAX ) THEN 
M = I 

LAMMAX = LAMHAT(I) 

END IF 
CONTINUE 



/* RATIOl IS FOR WITHOUT JACKNIFE */ 

/* RATI02 IS FOR WITH JACKNIFE */ 

DO 60 I = 1, N(II) 

RATIOl(I) = LAMBDA(I) / BIGLAM 
RATI02(I) = LAMHAT(I) / LAMMAX 
CONTINUE 



/* PART OF JACKNIFE METHOD FOR LAMBDA */ 

/* WITH OMIT 1 VARIABLE EACH TIME 

DO 90 I = 1, N(II) 

DO 80 J = 1, SAMPLE(I) 

LOMIT(I,J) = 0. 

DO 70 K = 1, SAMPLE(I) 

IF (J .NE. K) THEN 

LOMIT(I,J) = LOMIT(I,J) + CTIME(I,K) 
END IF 

CONTINUE 

LOMIT(I,J) = (SAMPLE(I)-l) / LOMIT(I,J) 

CONTINUE 

CONTINUE 
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/* ADAPT ABOVE RESULT TO OUR EQUATION TO */ 
/* GET THE RELIABILITY AND TRUE CONFIDENCE */ 
/* LIMIT */ 



SUMl = 0. 
SUM2 = 0. 
SUMS = 0. 
SUM4 = 0. 



DO 110 I = 1, N(II) 



/* NON JACKNIFING ( ORIGINAL ) */ 

SUMl = SUMl + RATIOl(I) * T(I) 

SUM2 = SUM2 + RATIOl(I) * TIME(I) 



/* WITH JACKNIFING */ 

ROMSUM = 0. 

DO 100 J = 1, NSTAR 

R0MIT(I,J) = LOMIT(I,J) / LOMIT(M,J) 
ROMSUM = ROMSUM + R0MIT(I,J) 

100 CONTINUE 



RSTAR(I) = NSTAR * RATI02(I) - 

* (NSTAR - 1) * ROMSUM / NSTAR 

SUMS = SUMS + RSTAR(I) * T(I) 

SUM4 = SUM4 + RSTAR(I) * TIME(I) 

110 CONTINUE 



/* R1 ; RELIABILITY OF A SYSTEM WITHOUT */ 
/* JACKNIFE */ 
/* R2 ; RELIABILITY OF A SYSTEM WITH */ 
/* JACKNIFE */ 



LAMHM = CHISQR / (2 * SUMl) 

R1(L) = EXP( - LAMHM * SUM2 ) 

LAMHST = CHISQR / (2 * SUMS) 
R2(L) = EXP( - LAMHST * SUM4 ) 



IF ( ABS(R1(L) - R2(D) . GE. DIFF(CASE) ) THEN 
DIFF(CASE) = ABS(RKL) - R2(L)) 
DIFFRl(CASE) = R1(L) 

DIFFR2(CASE) = R2(L) 

ENDIF 

120 CONTINUE 

/* NONIMSL LIBRARY 'SHSORT' WILL SORT Rl, R2 */ 
/* BY SHELL SORT ALGORITHM */ 
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CALL SHS0RT(R1, KEYl, NN) 
CALL SHS0RT(R2, KEY2, NN) 

MM = NN * (1 - ALPHA(LL)) 

RVALl(CASE) = Rl(MM) 
RVAL2(CASE) = R2(MM) 



/* SUBROUTINE FINDJ FINDS THE INDEX OF Rl, R2 */ 

/* WHICH THE VALUE OF IT IS CLOSEST TO RSYS. */ 

CALL FINDJ(R1, NN, RSYS(JJ), Jl) 

CALL FINDJ(R2, NN, RSYS(JJ), J2) 

CLRl(CASE) = Jl / FLOAT(NN) 

CLR2(CASE) = J2 / FLOAT(NN) 

130 CONTINUE 

140 CONTINUE 

150 CONTINUE 



WRITE(6,600) 

WRITE(6,650) 

WRITE(6,670) 

CASE = 1 

DO 210 II = 1, 2 
DO 200 JJ = 1, 2 
DO 190 LL = 1, 2 

WRITE(6,700) N(II), RSYS(JJ), ALPHA(LL), RVALl(CASE), 
* RVAL2(CASE), DIFF(CASE), CLRl(CASE), CLR2(CASE) 

WRITE(6,888) DIFFR1( CASE) , DIFFR2(CASE) 

CASE = CASE + 1 

190 CONTINUE 

200 CONTINUE 

WRITE(6,777) ( SAMPLE(J), J=1,N(II) ) 

WRITE(6,999) ( R(J), J=1,N(II) ) 

210 CONTINUE 



600 FORMATC ’ l' ,5(/) .7X, ***** RELIABILITY OF SERIES SYSTEM ****') 

650 FORMATC///, T50, 'ri, CLl ; WITHOUT JACKNIFING* , 

/,T50.'R2, CL2 ; WITH JACKNIFING') 

670 FORMATC ///,T6, ^NUMBER OF* ,T19 ,' RELIABILITY* ,T33, 'ALPHA* ,T46,*R1* , 

* T56 *R2' ,T63, 'MAX CRl - R2) * ,T79 ,* TRUE * ,T87 * TRUE* ,/ , 

* T6, COMPONENTS* ,T19, 'OF SYSTEM* ,T79, *C. L. 1\T87,*C.L. 2*, 

* /,T5,89C’-')) 



31 



o o o o 



700 FORMAT(/,T8,I5,T22,F5. 3,T33,F5. 2,T43,F7. 4,T53,F7. 4,T65,F8. 4,T78, 
* F7. 4,T86,F7. 4) 

777 FORMAT(/,T3, 'SAMPLE SIZE ARE 5(2X,I3) ) 

888 FORMAT(T60, 'Rl=' ,F6. 4, ' R2=',F6.4) 

999 FORMAT(/,T3,'R(I) ; ' ,5( 1X,F5. 3) ) 

STOP 

END 



* * 



» V V ^1^ ^v ^v ^ V ^ t 
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SUBROUTINE ZTABLE( ALPHA, ZALPHA) 

« SUBROUTINE ZTABL2 COMPUTES RIGHT PERCENT POINT ZALPHA FROM » 
« RIGHT CUMULATIVE PROBABILITY ALPHA » 

REAL ALPHA, ZALPHA 

IF ((ALPHA .GT. 0.0) .OR. (ALPHA . LT. 1)) THEN 
W = - AL0G(4 * ALPHA * ( 1 - ALPHA )) 

ZALPHA = SQRT(W * (2.06118 - ( 5.72622 / (W + 11.6406)))) 

IF (ALPHA .GT. 0.5) THEN 
ZALPHA = - ZALPHA 

ENDIF 

END IF 

RETURN 

END 






SUBROUTINE SHSORT( A,KEY,N) 
DIMENSION A(N),KEY(N) 

Ml=l 

6 M1=M1*2 

IF (Ml . LE. N) GO TO 6 

Ml=Ml/2-l 

MM=.MAX0( Ml/2,1) 

GO TO 21 

20 MM=MM/2 

IF (MM .LE. 0) GO TO 100 

21 K=N-MM 
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22 DO 1 J=1,K 
II=J 

11 IM=II+MM 

IF (A(IM) .GE. A(II)) GO TO 1 
TEMP=A(II) 

IT=KEY(II) 

A(II)=A(IM) 

KEY(II)=KEY(IM) 

A(IM)=TEMP 

KEY(IM)=IT 

IF (II .GT. 0) GO TO 11 
1 CONTINUE 
GO TO 20 
100 RETURN 
END 






VrVfVr^VVrVr'5VVrVrVc*7VVrVrVf‘V^VVr^V*Vr*^V^V^VVr*‘>V*^V^VV?Vf*VfVr/V'>V‘>V^V*'>VVrVr^VVr**VryfVf*’jV'>VV?VfV?VfV?'>VVcyf‘>VVr*VrVfVr^VVfVf'jV 



>> 

>> 



VALUE = ABS(A(NN) - R) 

DO 100 I = NN-1, 1, -1 

IF (ABS(A(I) - R) .LT. VALUE) THEN 
VALUE = ABS(A(I) - R) 

ELSE 

J = I + 1 
RETURN 

ENDIF 

100 CONTINUE 

RETURN 

END 



SUBROUTINE FINDJ(A, NN, R, J) 

« SUBROUTINE FINDJ FINDS THE INDEX OF ARRAY A WHICH THE 
« VALUE OF IT IS CLOSEST TO R. 

REAL A(NN), R, VALUE 
INTEGER J 
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APPENDIX C. FORTRAN CODE FOR INTERVAL ESTIMATION 
PROCEDURE . WEIBULL CASE 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



PROGRAM WEIBUL 



yrV?’>VVrVr'5VVf'5V':VV»'5VVc;^’‘V'5V'5Vyc'>V'5V‘V'5V'>VVc‘5V‘>VVc’'5V')V‘5V'>V'5VV? 



■jV'jVVc'VcVc'/VVc'jV'jV'jV'jV'jVVcVcVc’VrVfvV^^V^Vc’VoWc’jV'jWoV'jWf'sV^V 



* Vf 

* THIS PROGRAM COMPUTES THE RELIABILITY OF A SERIES SYSTEM * 

* USING TWO DIFFERENT METHODS FOR THE CASE WHEN FAILURE * 

* OF EACH COMPONENTS IS WE I BULLY DISTRIBUTED. * 

* AND ALSO COMPARES THE ESTIMATE OF SHAPE PARAMETERS ( BETA ) * 

* WHICH IS COMPUTED FROM TWO DIFFERENT METHODS. * 

■>v * 

* THE SERIES SYSTEM CONSIDERED IN THIS PORGRAM HAS * 

* N COMPONENT TYPES AND EACH COMPONENT HAS SAME NUMBER OF * 

* SAMPLE SIZES. * 



* BELOWS ARE GIVEN OR ASSUMED ; 

Vr * 

* BETA ; SHAPE PARAMETER OF WEIBULL DISTRIBUTION 

* ALPHA ; SIGNIFICANCE LEVEL 

* X ; TEST TIME * 

* RSYS ; RELIABILITY OF A SERIES SYSTEM * 

* 

* THESE ARE VARIABLES USED ; * 



* 



* 



* 






BETHAT 

BETNEW 

CHISQR 

CLR 

CLRNEW 

EXPO 

KSUBM 

LAMBDA 

LAMB I G 

LAMHAT 

LAMMAX 

LAMMU 

LAMSTR 

LAMNEW 

LAMUNR 

RATIO 

RATNEW 

RHAT 

RHTNEW 

XBENEW 

XHAT 

XIJHAT 

XIJNEW 

XNEW 

XTOBET 

W.WEIB 



BETA HAT FROM THE FORMULA 
BETA HAT FROM NEWTON - RAPHSON METHOD 
CHI-SQUARE VALUE FOR GIVEN ALPHA AND DF 
TRUE CONFIDENCE LIMIT FROM THE FORMULA 
TRUE CONFIDENCE LIMIT FROM THE N-R METHOD 
EXPONENTIAL RA.NDOM NUMBERS WITH MU=1 
K SUB M VALUE 

SCALE PARAMETER OF WEIBULL DISTRIBUTION 
MAXIMUM VALUE OF LAMBDA FROM THE N-R METHOD 
ESTIMATES OF LAMBDA FROM THE FORMULA 
MAXIMUM VALUE OF LAMBDA FROM THE FORMULA 
FINALLY COMPUTED LAMBDA FROM THE FORMULA 
LAMBDA TO THE BETA 

ESTIMATES OF LAMBDA FROM THE N-R METHOD 

FINALLY COMPUTED LAMBDA FROM THE N-R METHOD 

RATIO OF LAMBDA FOR THE FORMULA 

RATIO OF LAMBDA FOR THE N-R METHOD 

R( 500 * (1-ALPHA) ) FOR THE FORMULA 

R( 500 * (1-ALPHA) ) FOR THE N-R METHOD 

X TO THE BETNEW 

SUM OF XIJHAT 

W TO THE BETHAT 

W TO THE BETNEW 

SUM OF XIJNEW 

X TO THE BETHAT 

WEIBULL RANDOM NUMBERS 



* 



•iV 



y- 



* 

* 

Vr 
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ZALPHA ; RIGHT PERCENTILE POINT (NORMAL DISTRIBUTION) 



* 






PARAMETER ( N=5 , M=30 , NN= 500 ) 

PARAMETER ( N=15, M=30 , NN= 500) 

REAL BETA(N), LAMBDA(N), X(N), XHAT(N), XTOBET(N) , XX(N,M) 

REAL W(N,M), Y(N,M), BETHAT(N), WBHAT(N,M), LAMHAT(N) 

REAL RATIO(N), LAMSTR(N), LAMMAX, RSYS(2), KSUBM(N), LAMMU 

REAL EXPO(M), ALPHA(2), ZALPHA(2), RHAT(IO), RSL(NN), CHISQR 

REAL B(M), KEY(M), KEYl(NN), KEY2(NN), XIJHAT(N,M) 

REAL WEIB(M), BETNEW(N), XBENEW(N), XNEW(N), XIJNEW(N,M) 

REAL LAMNEW(N), LAMBIG, RATNEW(N), LAMUNR, RSLNEW(NN) 

REAL RHTNEW(IO), CLR(IO), CLRNEW(IO) 

REAL BHTBAR, BNRBAR, BHTMSE, BNRMSE 

INTEGER S(N), DF, ICOUNT, SAMPLE(N), CASE 

DATA ISEED / 1736 / 

DATA ALPHA / . 2 , .05 / 

DATA RSYS / . 90, . 975 / 

DATA BETA / 1. 2 , 1. 4 , 1. 6 , 1. 8 , 1. 5 / 

DATA BETA / 1. 2 , 1. 4 , 1. 6 , 1. 8 , 1. 5 , 1. 3 , 1. 7 , 1. 5 , 1. 2 , 1. 6 , 1. 9 , 1. 8 , 1. 4 , 

1. 6,1.1/ 

DATA BETA / 2. 2 , 2. 4 , 2. 6 , 2. 8 , 2. 5 / 

DATA BETA /2. 2 , 2. 4 , 2. 6 , 2. 8 , 2. 5 , 2. 3 , 2. 7 , 2. 5 , 2. 2 , 2. 6 , 2. 9 , 2. 8 , 2. 4 , 
2.6,2.!/ 

DATA X / 2,3,7,8,10/ 

DATA X /2, 3, 7, 8, 10, 3, 7, 10, 1,7, 8, 3, 10, 1,8/ 

DATA SAMPLE /9, 7, 10, 8, 6/ 

DATA SAMPLE /30 , 21 , 10 , 15 , 26/ 

DATA SAMPLE / 30,63,75,98,26 / 

DATA SAMPLE /6 , 7 ,5 , 6 , 7 , 8 ,9 ,5 ,8 , 7 ,9 , 10 , 7 , 9 , 6/ 

DATA SAMPLE /3 , 7 , 10 , 15 , 20 , 15 , 7 ,5 , 30 , 20 , 10 , 7 , 13 , 21 , 30/ 

DATA SA.MPLE / 15 ,40 , 35 , 17 , 26 , 67 , 50 , 65 , 80 , 32 , 95 , 100 , 15 ,45 , 30/ 

DATA SAMPLE / 15 * 15 / 

CASE = 1 



DO 10 I = 1, N 

IF ( SAMPLE(I) .LE. 15) KSUBM(I) = 1.40 
IF ( SAMPLE(I) .LE. 30) KSUBM(I) = 1.50 
10 CONTINUE 



DO 220 JJ = 1, 2 

C /•'^ FINDING THE LAMBDA AND LAMBDA STAR FROM THE GIVEN DATA */ 

DO 20 I = 1, N 
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LAMBDA(I) = (( - ALOG(RSYS(JJ)) / N) ** ( 1. /BETA( I) ) )/ X(I) 
LAMSTR(I) = LANBDA(I) ** BETA(I) 

S(I) =0.84 * SAMPLE(I) -0.5 

20 CONTINUE 

/* SUBROUTINE ZTABLE COMPUTES THE RIGHT PERCENT POINT ZALPHA */ 
/* FROM RIGHT CUMULATIVE PROBABILITY ALPHA */ 

DO 40 I = 1, 2 

CALL ZTABLE(ALPHAd) , ZALPHA(I)) 

40 CONTINUE 



C /* COMPUTE THE DEGREE OF FREEDOM FOR CHI-SQUARE */ 

DF = 0 

DO 60 I = 1, N 

DF = DF + SAMPLE(I) 

60 CONTINUE 

DF = 2 * DF 



C /* COMPUTE CHI-SQUARE VALUE FOR GIVEN ALPHA AND DF */ 

DO 200 LL = 1, 2 

IF ( DF .EQ. 1 ) THEN 
P = ALPHA(LL) / 2 
CALL ZTABLECP, ZALPHA(LL)) 

CHISQR = ZALPHA(LL) ** 2 

ELSE IF ( DF .EQ. 2 ) THEN 

CHISQR = -2 * ALOGCALPHA(LL)) 

ELSE IF ( DF .GE. 3 ) THEN 

CHISQR = DF * (1 - 2. /(9 * DF) + 

* ZALPHA(LL) * SQRT( 2./(9 * DF))) ** 3 

END IF 

LAMMAX = 0 

LAMB I G = 0 

DO 180 ITER = 1, NN 
DO 140 I = 1, N 

/* GENERATE EXPONENTIAL RANDOM NUMBER WITH MU=1 */ 
CALL LEXPNCISEED, EXPO, SAMPLE(I), 1, 0) 

/* CONVERT EXPONENTIAL RANDOM NUMBER WITH MU=1 */ 
/* TO MU=LAMBDA AND GET WEIBULL RANDOM NUMBERS */ 
DO 80 J = 1, SAMPLE(I) 

XX(I,J) = EXPO(J) / LAMSTR(I) 

W(I,J) = XX(I,J) ** (1/BETA(I)) 

Y(I,J) = ALOG(W(I,J)) 

WEIB(J) = W(I,J) 
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B(J) = Y(I,J) 

80 CONTINUE 

/* SUBROUTINE SHSORT WILL SORT ARRAY B IN 
/* ASCENDING ORDER BY SHELL SORT ALGORITHM 
CALL SHSORTCB, KEY, SAMPLE(I)) 

/* FINDING BETA HAT (BETHAT) BY THE FORMULA 

SUMl = 0 
SUM2 = 0 

DO 100 J = 1, S(I) 

SUMl = SUMl + B(J) 

100 CONTINUE 

DO 120 J = S(I)+1, SAMPLE(I) 

SUM2 = SUM2 + B(J) 

120 CONTINUE 

BETHAT(I) = (SAMPLE(D) * KSUBM(I) / 

* (S(I)*SUM2 / (SAMPLE(I) - S(I)) - SUMl) 

TEMP = BETA(I) 

/* SUBROUTINE NEWTON WILL COMPUTE THE ESTIMATE 
/* OF BETA (BETNEW) BY NEWTON-RAPHSON METHOD 
CALL NEWTONCWEIB, SAMPLE(I), TEMP, BETNEW(I)) 

XTOBET(I) = X(I) ** BETHAT(I) 

XBENEW(I) = X(I) ** BETNEW(I) 

XHAT(I) = 0. 

XNEW(I) = 0. 

DO 130 J = 1, SAMPLE(I) 

XIJHAT(I,J) = W(I,J) ** BETHAT(I) 

XHAT(I) = XHAT(I) + XIJHAT(I,J) 

XIJNEW(I,J) = W(I,J) ** BETNEW(I) 

XNEW(I) = XNEW(I) + XIJNEW(I,J) 

130 CONTINUE 

LAMHAT(I) = SAMPLE(I) / XHAT(I) 

LAMNEW(I) = SAMPLE(I) / XNEW(I) 

IF (LAMHAT(I) . GT. LAMMAX ) THEN 
LAMMAX = LAMHAT(I) 

END IF 

IF (LAMNEW(I) .GT. LAMB I G ) THEN 
LAMBIG = LAMNEW(I) 

END IF 

140 CONTINUE 

SUM3 = 0 
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SUM3N = 0 
SUM4 = 0 
SUM4N = 0 

DO 160 I = 1, N 

RATIO(I) = LAMHAT(I) / LAMMAX 
RATNEW(I) = LAMNEW(I) / LAMB I G 
SUM3 = SUM3 + RATIO(I) * XHAT(I) 

SUM4 = SUM4 + RATIO(I) * XTOBET(I) 

SUM3N = SUM3N + RATNEW(I) * XNEW(I) 
SUM4N = SUM4N + RATNEW(I) * XBENEW(I) 

160 CONTINUE 

LAMMU = CHISQR / ( 2 * SUM3 ) 

RSL(ITER) = EXP( - LAMMU * SUM4 ) 

LAMUNR = CHISQR / ( 2 * SUM3N ) 

RSLNEW(ITER) = EXP( - LAMUNR * SUM4N ) 

180 CONTINUE 

CALL SHSORT(RSL, KEYl, NN) 

CALL SHSORT(RSLNEW, KEY2, NN) 

KK = ( 1 - ALPHA(LL)) * NN 



/* RHAT IS THE RELIABILITY OF THE SERIES SYSTEM COMPUTED */ 
/* BY THE FORMULA */ 

/•■ RHTNEV IS THE RELIABILITY OF THE SERIES SYSTEM COMPUTED */ 
/* BY THE NEWTON - RAPHSON METHOD */ 

RHAT(CASE) = RSL(KK) 

RHTNEW(CASE) = RSLNEW(KK) 

/•'• SUBROUTINE FINDJ WILL FIND THE TRUE CONFIDENCE LIMIT */ 
CALL FINDJ(RSL, NN, RSYS(JJ), Jl) 

CALL FINDJ(RSLNEW, NN, RSYS(JJ), J2) 



CLR(CASE) = Jl / FLOAT(NN) 

CLRNEW(CASE) = J2 / FLOAT(NN) 

CASE = CASE + 1 



200 CONTINUE 



220 CONTINUE 



DO 240 I = 1, N 

WRITE(6,500) BETA(I), SAMPLE(I) 
240 CONTINUE 



WRITE(6,550) N, NN 
WRITE(6,600) 
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WRITE(6,650) 

WRITE(6,670) 



CASE = 1 

DO 280 JJ = 1, 2 

DO 260 LL = 1, 2 

WRITE(6,700) N,RSYS(JJ) ,ALPHA(LL) , RHAT( CASE ) ,RHTNEW( CASE) , 
* CLR(CASE), CLRNEW(CASE) 

CASE = CASE + 1 

260 CONTINUE 

280 CONTINUE 



500 FORMATCTIO, ' INITIAL BETA :’,3X,F4. 1, 

* T45, 'NUMBER OF SAMPLE SIZES ;',2X,I3) 

550 FORMATCTIO, 'NUMBER OF COMPONENTS ;',2X,I3, 

* /,T10, 'NUMBER OF REPLICATIONS :',2X,I3) 

600 FORMATC ' 1' ,5(/),7X, '**** RELIABILITY OF SERIES SYSTEM ****') 

650 FORMATC ///, T35, 'r 1, CLl ; USING THE FORMULAE', 

/,T35 'R2, CL2 ; USING THE NEWTON RAHPSON METHOD') 

670 FORMATC ///,T6, NUMBER OF' T19 , 'RELIABILITY' ,T33, 'ALPHA' ,T46, 'Rl' , 

* T56 'R2' ,T66, 'true' ,T74, 'TRUE' ,/, 

* T6. 'components' ,T19, 'OF system', T 66,'C.L. l',T74,'C.L. 2', 
/,T5,76C'-')) 

700 FORMATC/,T8,I5,T22,F5. 3,T33,F5. 2,T43,F7. 4,T53,F7. 4,T65, 

* F7. 4,T73,F7. 4) 



STOP 

END 




SUBROUTINE ZTABLEC ALPHA , ZALPHA) 

« SUBROUTINE ZTABL2 COMPUTES RIGHT PERCENT POINT ZALPHA FROM » 
« RIGHT CUMULATIVE PROBABILITY ALPHA » 

REAL ALPHA, ZALPHA 



IF CCALPHA .GT. 0.0) .OR. CALPHA . LT. 1)) THEN 
W = - AL0GC4 * ALPHA * C 1 - ALPHA )) 

ZALPHA = SQRTCW * C 2. 06118 - C 5.72622 / CW + 11.6406)))) 

IF CALPHA .GT. 0.5) THEN 
ZALPHA = - ZALPHA 

END IF 
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ENDIF 



RETURN 

END 
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SUBROUTINE SHSORT(A,KEY,N) 
DIMENSION A(N),KEY(N) 

Ml=l 

6 M1=M1*2 

IF (Ml .LE. N) GO TO 6 
Ml=Ml/2-l 
MM=MAX0(Ml/2, 1) 

GO TO 21 

20 MM=MM/2 

IF (MM .LE. 0) GO TO 100 

21 K=N-MM 

22 DO 1 J=1,K 
II=J 

11 IM=II+MM 

IF (A(IM) .GE. A(II)) GO TO 1 
TEMP=A(II) 

IT=KEY(II) 

A(II)=A(IM) 

KEY(II)=KEY(IM) 

A(IM)=TEMP 

KEY(IM)=IT 

II=II-MM 

IF (II . GT. 0) GO TO 11 
1 CONTINUE 
GO TO 20 
100 RETURN 
END 



* 



SUBROUTINE NEWTON(TIME, M, BETA, BETNEW ) 

« SUBROUTINE NEWTON WILL COMPUTE THE ESTIMATE OF BETA (BETHAT) » 
« BY NEWTON - RAPHSON METHOD » 

REAL TIME(M), LNTIME(IOO), A(3), LNTS 

INTEGER R 

C DATA TS / 6 / 

C =0 
ITER = 0 
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R = M 

C LNTS = ALOG(TS) 

DO 20 I = 1, R 

LNTIME(I) = ALOG(TIMEd)) 
C = C + LNTIME(I) 

20 CONTINUE 



C = C / R 

30 DO 60 J = 1, 3 
SUM = 0 



DO 40 K = 1, R 

IF ( J .EQ. 1 ) THEN 

SUM = SUM + TIME(K) ** BETA 

ELSE 

SUM = SUM + ((TIME(K) ** BETA)*(LNTIME(K) ** (J-1))) 
ENDIF 

40 CONTINUE 



A(J) = SUM 
60 CONTINUE 



/* FUNCTION FPRIME IS THE DERIVATIVES OF FUNCTION F */ 

QUOT = A(2) / A(l) 

FPRIME = A(3) / A(l) - (QUOT**2) + ( ( 1/BETA)**2) 

F = QUOT - (1 / BETA) - C 

BETA IS UPDATED EACH TIME AND CHECK IF IT CONVERGES */ 

BETA = BETA - F / FPRIME 

ITER = ITER + 1 

IF (BETA .GT. 25.) GOTO 100 

IF ( ABS(F) .GT. 0.0001 ) GOTO 30 

ALPHA = (A(l) / R) ** (1/BETA) 

BETNEW = BETA 



RETURN 

100 VRITE(6,*) 'DID NOT CON\TRGE' 

WRITE(6,*) 'TRY AGAIN WITH BETTER ESTIMATE OF BETA' 

RETURN 

END 









^10 , 
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SUBROUTINE FINDJ(A, NN, R, J) 

« SUBROUTINE FINDJ FINDS THE INDEX OF ARRAY A WHICH THE » 

« VALUE OF IT IS CLOSEST TO R. » 

REAL A(NN), R, VALUE 
INTEGER J 

VALUE = ABS(A(NN) - R) 

DO 100 I = NN-1, 1, -1 

IF (ABS(A(I) - R) .LT. VALUE) THEN 
VALUE = ABS(A(I) - R) 

ELSE 

J = I + 1 
RETURN 

END IF 

100 CONTINUE 

RETURN 

END 
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APPENDIX D. FORTRAN CODE FOR INTERVAL ESTIMATION 
PROCEDURE - NORMAL CASE 



PROGRAM NORMAL 






* THIS PROGRAM DETERMINE THE ACCURACY OF AN APPROXIMATE 

* LOWER CONFIDENCE BOUND FOR P( X > Y ). * 

“V Vr 

* PROGRAM NORMAL IS THE CASE WHEN Y IS GIVEN A VALUE YO * 

* & X IS NORMALLY DISTRIBUTED WITH * 

UNKNOWN PARAMETERS. * 

* Vf 



'.-***'!V***?V**'5V-V'>V***'>V*********Vf')V*****'5V****?V**->V*-V*****5Wr)V*-.V***>V*'5V**** 



REAL YO, P, ZP, MUX, SUMX, SUMX2 , XBAR 

REAL TEMPI, TEMP2, TEMP 3, SI GHAT 

REAL R(3), SIGMAX(3), ALPHA(2), X(IOO), Y(IOO) 

REAL RL1(54,1000), ARLl(lOOO), BRL1(54), KEY(IOOO), CIRL1(54) 
REAL T1RLK54), T1CI1(54), T2RL1(54), T2CI1(54) 

REAL TRLK2), TCI 1(2) 

REAL ZHAT(IOOO), ZKNIFE(54), ZVAR(54), ZBAR(54) 

REAL MUJUNK(54) 

INTEGER NUMX(3), NUMY(3), CASE, LINE 



DATA R /. 95 ,. 99 ,. 995/, SIGMAX /. 5 , 1 , 20/ , NUMX /10,25,75/ 

DATA NUMY /10,25,75/ , ALPHA / .2, .1/ 

DATA ISEED, JSEED / 4875,7981 / 

DATA YO / 400 / 

CASE = 0 

/* II IS THE INDEX FOR R. R(l) = .95, R(2) = .99, R(3) = .995 */ 

DO 500 II = 1, 3 



/* P IS THE RIGTH(UPPER) CUMULATIVE PROBABILITY */ 

P = 1 - R(II) 

/* SUBROUTINE ZTABL2 COMPUTES RIGHT PERCENT POINT ZP */ 

/* FROM RIGHT CUMULATIVE PROBABILITY P */ 

CALL ZTABL2C P, ZP ) 

/* JJ IS THE INT3EX FOR SIGMA OF X. 

/* SIGMAX(l) = .5, SIGMAX(2) = 1, SIGMAX(3) = 3 */ 
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DO 400 JJ = 1, 3 

MUX = ZP * SIGMAX(JJ) + YO 



/* KK IS THE INDEX OF NUMBER OF X. */ 
/* NUMX(l) = 10, NUMX(2) = 25, NUMX(3) = 75 */ 
DO 300 KK = 1, 3 



/* LL IS THE INDEX FOR ALPHA. */ 

/* ALPHA(l) = .2, ALPHA(2) = .1 */ 

DO 200 LL = 1, 2 

CALL ZTABL2(ALPHA(LL) , ZALPHA) 

rA<?F = rA9P + 1 

TEMP = (NUMX(KK) - 1. ) / NUMX(KK) 

TEMPO = SQRT(TEMP) 

MUJUNK(CASE) = MUX 



/* REPLICATE 1000 TIMES FOR EACH CASES */ 

DO 100 I = 1, 1000 

/* USE NORMAL RANDOM NUMBER GENERATOR TO GET */ 

/* NUMX(KK) NUMBER OF X. */ 

CALL LN0RM( I SEED , X , NUMX( KK) , 2 , 0 ) 



/* THIS PART IS TO GET SAMPLE MEAN(XBAR) AND */ 

SAMPLE VARIANCE(XVAR) */ 



DO 50 MM = 1, NUMX(KK) 

X(MM) = SIGMAX(JJ) * X(MM) + MUX 
50 CONTINUE 



/* SUBROUTINE VAR WILL COMPUTE SAMPLE MEAN AND */ 

/* SAMPLE VARIANCE */ 

CALL VAR(X, NUMX(KK), XBAR, XVAR) 

/* NOW WE COMPUTE THE LOWER CONFIDENCE BOUND */ 

TEMPI = (XBAR - YO) / SQRT(XVAR * TEMPO) 

TEMP2 = 1. / NUMX(KK) + 

((XBAR - YO) ** 2) / (2*(NUMX(KK)+1)* XVAR) 

TEMPI = TEMPI - ZALPHA * SQRT(TEMP2) * TEMP 



/* SUBROUTINE ZTABLl COMPUTE RIGHT CUMULATIVE PROBA- */ 
/* BILITY FROM RIGHT PERCENT POINT. */ 

CALL ZTABLKTEMPl, ARLl(I)) 
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200 

300 

400 

500 



ARLl(I) = 1. - ARLl(I) 

CALL JKNIFECX, NUMX(KK), YO, ZHAT(I)) 

CONTINUE 

CALL VARCZHAT, 1000, ZBAR(CASE), ZVAR(CASE)) 

SIGHAT = l./NUMX(KK) + ZBAR(CASE)**2 / ( 2*(NUMX(KK)+1) ) 

SIGHAT = SQRT( SIGHAT) * TEMP 

SIGHAT = ZBAR(CASE) / TEMPO - SIGHAT * ZALPHA 

CALL ZTABLK SIGHAT, ZKNIFE(CASE) ) 

ZKNIFE(CASE) = 1 - ZKNIFE(CASE) 

/* NON-IMSL LIBRARY 'SHSORT' WILL SORT ARLl BY SHELL 
/* SORT ALGORITHM. 

CALL SHS0RT(ARL1, KEY, 1000) 

DO 150 I = 1, 1000 

RL1(CASE,I) = ARLl(I) 

CONTINUE 



/* SUBROUTINE FINDJ FINDS THE INDEX OF ARL20 WHICH THE 
/* VALUE OF IT IS CLOSEST TO R. 

CALL FINDJ(ARL1, R(II), J) 



/* DIVIDING THE INDEX BY 1000 WILL GIVE US TRUE 
/* CONFIDENCE LEVEL. 

CIRLl(CASE) = J / 1000. 

MM = 1000 * (1 - ALPHA(LL)) 

BRLl(CASE) = RL1(CASE,MM) 

CALL TNYVAL(R(II) , SIGMAX(JJ), NUMX(KK), 

ALPHA(LL), ZALPHA, MUX, YO, TRLl, TCIl) 

TlRLl(CASE) = TRLK 1) 

TlCIl(CASE) = TCIl(l) 

T2RLKCASE) = TRL1(2) 

T2CIKCASE) = TCIK2) 



CO.NTINUE 

CO.NTINUE 

CONTINUE 

CONTINUE 



WRITE(6,600) 

WRITE(6,630) 

I = 1 
LINE = 0 

DO 550 II = 1, 3 
DO 550 JJ = 1, 3 
DO 550 KK = 1, 3 

IF (LINE .GE. 3) THEN 
WRITE(6,800) 

LINE = 0 
ENDIF 



WRITE(6,710) R(II), SIGMAX(JJ), NUMX(KK), 

* BRLl(I), BRLKI+I), CIRLl(I), CIRL1(I+1), 

* ZKNIFE(I), ZKNIFE(I+1), 

* MUJUNKC I ) , MUJUNKC I+l) 

1 = 1 + 2 

LINE = LINE + 1 
550 CONTINUE 



WRITE(6,610) 

WRITE(6,650) 

I = 1 
LINE = 0 

DO 560 II = 1, 3 
DO 560 JJ = 1, 3 
DO 560 KK = 1, 3 

IF (LINE .GE. 3) THEN 
WRITE(6,800) 

LINE = 0 
ENDIF 

WRITE(6,700) R(II), SIGMAX(JJ), NUMX(KK), 

* TlRLl(I), TlRLKI+1), TlCIl(I), T1CI1(I+1) 

1 = 1 + 2 
LINE = LINE + 1 
560 CONTINUE 



WRITE (6, 620) 

WRITE(6,650) 

I = 1 
LINE = 0 

DO 570 II = 1, 3 
DO 570 JJ = 1, 3 
DO 570 KK = 1, 3 

IF (LINE .GE. 3) THEN 
WRITE (6, 800) 

LINE = 0 
ENDIF 
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WRITE(6,700) R(II), SIGMAX(JJ), NUMX(KK), 

* T2RLKI), T2RLKI+1), T2CI1(I), T2CI1(I + 1) 

1 = 1 + 2 
LIKE = LINE + 1 
570 CONTINUE 



600 FORMATC 1' ,5(/),7X,'**** Y IS GIVEN A VALUE OF YO ****') 

610 FORMATC 1' ,5(/),7X,’**** Y IS GIVEN A VALUE OF YO ****', 

* //,20X,'-- X IS TRUNCATED NORMAL WITH PROBABILITY .95 
620 FORMATC’ l',5(/),7X,’**** Y IS GIVEN A VALUE OF YO ****’ ^ 

* //,20X,'-- X IS TRUNCATED NORMAL WITH PROBABILITY .90 

630 FORMATC 3(/),5X,’ CASE ',6X,' RLlC 1000*C 1-ALPHA) ) 

* 4X,’TRUE CONFIDENCE LEVEL’ ,5X, 'JACKNIFE' , 

///,5X,’R SX N’ ,3X,’ALPHA = . 2’ ,12X,’. 1’ ,11X,’. 2' , 

* 12X,’. 1’ ,/,3X,13C'-’ ),6X,25C’-’) CX,2lC-'),2X,16C ’-'),/) 

650 FORMATC3C/) ,5X, ’ CASE ’,6X,’ RLlC 1000*C 1 -ALPHA) ) ’, 

* 4X, ’TRUE CONFIDENCE LEVEL’ ,///,5X, ’R SX N’,3X. 

* ’ALPHA = .2’ ,12X.’. 1’ ,11X,’.2’ ,12X,’. 1’ ,/,3X,13C'- ),6X, 

25C'-’),4X,21C’-C,/) 

700 FORMATC 3X,F4. 3,2X,F3. 1 , 2X , 12 , 4C 7X ,F7. 4) ) 

710 FORMATC 3X,F4. 3, IX, F4. 1 , 2X , 12 ,4C 7X,F7. 4) , IX, 2C 2X , F6. 4) ,2C2X,F9. 5)) 
800 FORMATC/) 

RETURN 

END 
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SUBROUTINE ZTABLlCZALPHA, ALPHA) 

« SUBROUTINE ZTABLl COMPUTE RIGHT CUMULATIVE PROBABILITY FROM » 
« RIGHT PERCENT POINT. » 

REAL PI, ALPHA, Z ALPHA 

PARA.METER C PI = 3. 141592 ) 



IF C ZALPHA . EQ. 0. 0 ) THEN 
ALPHA = 0.5 

ELSE 

ALPHA = EXPC -2 * ZALPHA ** 2 / PI ) 

ALPHA = SQRTCl - ALPHA * Cl + 2*CPI-3)*ZALPHA**4/C3*PI**2)) ) 
ALPHA = 0.5 * C 1 - ALPHA ) 

IF C ZALPHA . LT. 0.0 ) THEN 
ALPHA = 1 - ALPHA 
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ENDIF 

ENDIF 

RETURN 

END 
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SUBROUTINE ZTABL2( ALPHA, ZALPHA) 

« SUBROUTINE ZTABL2 COMPUTES RIGHT PERCENT POINT ZALPHA FROM » 
« RIGHT CUMULATIVE PROBABILITY ALPHA » 

REAL ALPHA, ZALPHA 

IF ((ALPHA .GT. 0.0) .OR. (ALPHA . LT. 1)) THEN 
W = - AL0G(4 * ALPHA * ( 1 - ALPHA )) 

ZALPHA = SQRT(V * (2.06118 - ( 5.72622 / (V + 11.6406)))) 

IF (ALPHA .GT. 0.5) THEN 
ZALPHA = - ZALPHA 

ENDIF 

ENDIF 

RETURN 

E.ND 






SUBROUTINE TTABLE( ALPHA, NU, TALPHA) 

« SUBROUTINE TTABLE COMPUTES RIGHT CUMULATIVE PROBABILITY » 

« TALPHA FROM T-DISTRIBUTION BY GIVEN ALPHA AND NU » 

REAL ALPHA, ZALPHA, TALPHA, Al, A2, A3 
I.NTEGER NU 



IF ((ALPHA .GT. 0.0) .OR. (ALPHA . LT. 1)) THEN 
CALL ZTABL2( ALPHA, ZALPHA) 

Al = (ZALPHA**2 + 1) / 4 

A2 = (5*ZALPHA**4 + 16*ZALPHA**2 + 3) / 96 

A3 = (3*ZALPHA**6 + 19*ZALPHA**4 + 17*ZALPHA**2 - 15) / 384 

TALPHA = ZALPHA * (1 + Al/NU + A2/NU**2 + A3/NU**3) 



IF ( ALPHA .GT. 0.5 ) THEN 
TALPHA = - TALPHA 

ENDIF 



4S 
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ENDIF 



RETURN 

END 
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SUBROUTINE SHSORT( A,KEY,N) 
DIMENSION A(N),KEY(N) 

Ml=l 

6 M1=M1*2 

IF (Ml .LE. N) GO TO 6 

Ml=Ml/2-l 

MM=MAX0(M1/2,1) 

GO TO 21 

20 MM=MM/2 

IF (MM .LE. 0) GO TO 100 

21 K=N-MM 

22 DO 1 J=1,K 
II=J 

11 IM=II+MM 

IF (A(IM) . GE. A(II)) GO TO 1 
TEMP=A(II) 

IT=KEY(II) 

A(II)=A(IM) 

KEY(II)=KEY(IM) 

A(IM)=TEMP 

KEY(IM)=IT 

II=II-MM 

IF (II .GT. 0) GO TO 11 
1 CONTINUE 
GO TO 20 
100 RETURN 
END 



'kk' 









SUBROUTINE VAR(Z, NUMZ, ZBAR, ZVAR) 

« SUBROUTINE VAR WILL COMPUTE THE MEAN AND THE VARIANCE » 

« OF ARRAY Z » 

REAL Z(NUMZ), SUMZ, SUMZ2, ZBAR, ZVAR 

SUMZ = 0 
SUMZ2 = 0 

DO 100 I = 1, NUMZ 

SUMZ = SUMZ + Z(I) 
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SUMZ2 = SUMZ2 + Z(I) ** 2 



100 CONTINUE 

ZBAR = SUMZ / NUMZ 

ZVAR = ABS( SUMZ2 / NUMZ - ZBAR ** 2 ) 

RETURN 

END 
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SUBROUTINE FINDJ(A, R, J) 

« SUBROUTINE FINDJ FINDS THE INDEX OF ARRAY A WHICH THE » 

« VALUE OF IT IS CLOSEST TO R. » 

REAL A(IOOO), R, VALUE 
INTEGER J 

VALUE = ABS(AdOOO) - R) 

DO 100 I = 999, 1, -1 

IF (ABS(A(I) - R) .LT. VALUE) THEN 
VALUE = ABS(A(I) - R) 

ELSE 

J = I + 1 
RETURN 

ENDIF 

100 CONTINUE 

RET’ RN 
END 
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SUBROUTINE TNYVAL( R , S IGMAX , NUMX , ALPHA , ZALPHA , MUX , YO , 

TBRL1,TCIR1) 

« SUBROUTINE TNYVAL IS THE CASE WHEN X IS TRUNCATED NORMAL » 
« AND Y IS GIVEN A VALUE » 



REAL R, S IGMAX, ALPHA, TALPHA, MUX, P(2), ZP(2), A(2) 

REAL X(IOO), XI, TEMPI, TEMP2, YO 

REAL TARLl(lOOO), KEY(IOOO), TRLl(lOOO), TBRL1(2), TCIR1(2) 
INTEGER NUMX, CASE, COUNT 
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DATA P / .95, .90 / 

DATA ISEED, JSEED / 4875, 7981 / 

DO 400 I = 1, 2 

P(I) = 1 - P(I) 

COUNT = 0 

CALL ZTABL2(P(I), ZP(I)) 

A(I) = MUX + ZP(I) * SIGMAX 
DO 200 J = 1, 1000 

100 CALL LNORM( ISEED, XI, 1, 2, 0) 

XTEMP = MUX + SIGMAX * XI 

IF (XTEMP .LE. A(D) THEN 
COUNT = COUNT + 1 
X( COUNT) = XTEMP 

ELSE 

GOTO 100 
END IF 

IF (COUNT .LT. NUMX) GOTO 100 
CALL VAR(X, NUMX, XBAR, XVAR) 



/* NOW WE COMPUTE THE LOWER CONFIDENCE BOUND 

TEMPI = (XBAR - YO) / SQRT( XVAR* ( NUMX- 1 ) / NUMX) 
TEMP2 = 1. / NUMX + 

((XBAR - YO) ** 2) / (2*(NUMX+1) * XVAR) 
TEMPI = TEMPI - ZALPHA * SQRT(TEMP2) 



/* SUBROUTINE ZTABLl COMPUTE RIGHT CUMULATIVE PROBA- 
/* BILITY FROM RIGHT PERCENT POINT. 

CALL ZTABLKTEMPl, TARLl(J)) 

TARLl(J) = 1. - TARLl(J) 

200 CONTINUE 



/* NON-IMSL LIBRARY ’SHSORT’ WILL SORT TARLl BY SHELL 
/* SORT ALGORITHM. 

CALL SHSORT( TARLl, KEY, 1000) 

DO 300 J = 1, 1000 
TRLl(J) = TARLl(J) 



o n n o o o no 



300 CONTINUE 



/* SUBROUTINE FINDJ FINDS THE INDEX OF TARLl WHICH THE */ 

/* VALUE OF IT IS CLOSEST TO R. */ 

CALL FINDJCTARLl, R, J) 



/* DIVIDING THE INDEX BY 1000 WILL GIVE US TRUE */ 

/* CONFIDENCE LEVEL. */ 

TCIRl(I) = J / 1000. 

MM = 1000 * (1 - ALPHA) 

TBRLl(I) = TRLl(MM) 

400 CONTINUE 



RETURN 

END 
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SUBROUTINE JKNIFE(X, N, YO, ZHAT) 

« THIS ROUTINE IS FOR THE 'JACKNIFE' METHOD. » 

« DUMMY PARAMETER IS ZHAT. » 

REAL X(N) ,XOMIT( 100) ,XOMBAR( 100) ,XOMVAR( 100) , ZOMHAT( 100) ,ZHAT 
INTEGER M, N 

M = N - 1 
SUM = 0. 



DO 200 I = 1, N 

DO 100 J = 1, M 

IF (J .GE. I) THEN 

XOMIT(J) = X(J+1) 

ELSE 

XOMIT(J) = X(J) 

END IF 

100 CONTINUE 

CALL VAR(XOMIT, M, XOMBAR(I), XOMVAR(I)) 

ZOMHAT(I) = ( XOMBAR(I) - YO ) / SQRT( XOMVAR( I ) ) 
SUM = SUM +ZOMHAT(I) 



200 CONTINUE 



ZHAT = SUM / N 
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RETURN 

END 
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